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Abstract 

We develop some calculation schemes to determine dynamics of a wide class 
of integrable quantum-optical models using their symmetry adapted reformulation 
in terms of polynomial Lie algebras su p d(2). These schemes, based on "diagonal" 
representations of model evolution operators (via diagonalizing Hamiltonians with 
the help of the su p d(2) defining relations), are implemented in the form adapted 
for numerical calculations. Their efficiency is demonstrated on the example of the 
second-harmonic-generation model. 

1 Introduction 

During last decades a great attention is being paid to examine different quantum-optical 
models with Hamiltonians given by nonlinear functions in Lie algebra generators since 
they enable to reveal new physical effects and phenomena (see, e.g., [1-10] and references 
therein). To analyze such models one uses mainly numerical calculations [10] because 
standard Lie-algebraic techniques well adapted for solving problems with linear (in Lie 
algebra generators) Hamiltonians [11], are non-efficient here, and most of other analytical 
techniques (e.g., the algebraic Bethe ansatz [12]) require in general cases tedious calcu- 
lations and do not yield simple analytical expressions for physical quantities. On other 
hand, standard numerical calculation schemes dealing with initial formulations of models 
[10] are limited by computer powers and are not adapted to reveal many peculiarities of 
model dynamics [8,9]. 

However, recently a new universal Lie-algebraic approach, essentially improving both 
analytical and numerical solutions of physical problems, has been suggested in [4,5] and 
developed in [5-9] for some nonlinear quantum models whose Hamiltonians H have in- 
variance groups Gi : [Gi, H] = 0. It is based on reformulations of models under study 
in terms of (introduced in [4,5,13]) polynomial Lie algebras (PLA) g v ^ as dynamic sym- 
metry algebras g D : g D = g p d completely describing model dynamics. Specifically, two 
analytical approximations of evolution operators Unit) were found in [8] for a wide class 
of quantum-optical models with g D = g p( i = su p d(2). Furthermore, appropriate Lie- 
algebraic path- integral schemes for solving physical problems were developed in [9] . They 
enable to examine model dynamics at quasiclassical levels described by "large" values 
of certain characteristic parameters [8]. However, up to now systematic examinations 
of their accuracy and efficiency were absent because SM p d(2)-techniques of getting exact 
solutions [5-7] yield only difference and difference-differential equations and fairly com- 
plicated (for practical calculations) algorithms rather than closed analytical expressions 
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required for this aim. In the present work we cancel in part this lacune. Specifically, 
after some preliminaries (Section 2) we develop an algorithm and routines to implement 
an exact calculation scheme for determining diagonal representations of model evolution 
operators via solving a spectral problem by means of the s-u p( i(2)-techniques [5] (Section 
3). Their efficiency is demonstrated with the help of the computer experiments for the 
second-harmonic-generation model which are also used to investigate the accuracy of the 
appropriate quasiclassical approximation obtained in [8] (Section 4). 

2 Models and their symmetry adapted su p d(2) - for- 
mulations 

We consider a wide-spread class of quantum-optical models with Hamiltonians [1,3] 



H mps (m;n) = h 



Uiafai + ^0^0 °o + 9 a t ■ ■ ■ <4 n a o + 9* o>h ■ ■ ■ a in a o 



i=l 



1 < i\ < ?2 • • • < i n < m , n > 2 (2.1a) 



M 



H D = h^j 1 aia 1 + Y,[^/ 2 + 9<J+(i)(aiT + 9^-matT]^ (2.16) 

where operators of field modes with frequencies u h a a (i) are Pauli matri- 

ces, e is an energy difference of two-level atoms, g are coupling constants, and non- 
quadratic parts of Eqs. (2.1a) describe different processes of multiphoton scattering in- 
cluding non-degenerated frequency conversion (H mps (n; n)) and higher harmonics genera- 
tion (H mps (l;n)) while Eq. (2.1b) yields n-photon point-like Dicke models of matter- 
radiation interactions [1,8]. Hilbert spaces L(H mps ) are Fock spaces L F (k < m) — 
Span{\{m}) oc nf = oK + ) ni |0>} and L(H D ) = L F (1) ® L a , L a = Span^^ |±>(i)} where 
|±)(i)} are eigenstates of i-th atom [3,4]. 
According to [5,9] Hamiltonians (2.1) have invariance groups 

Gi = C n ®l[U j (l), C n = {a+ - exp(i27rA;/n)a+}, W(l) = {exp(^} (2.2) 

3 

where for H mps (m;n)Rj G Span{afai = Ni} and for H D Rj are functions in J^fti^aii) 
and, besides, Gi has the permutation group Sm = {& a (i) ~^ a a(j)} as an extra factor. 
Therefore, one can introduce (via the Gj- invariant polynomial Jordan mappings [5,9]) 
two sets of collective operators: integrals of motion Rj = Rj(af,di, . . .),j = 1, . . . and 
Gj-mvariant dynamic variables V a= o,± = V a= o,±(ai, a~j~, . . .) obeying the commutation re- 
lations 

Wo, V±] = ±V±, W-, V+] = <f>(V ; {Rj}) = 4>(V + 1; {Rj}) - tp(V ; {Rj}), 



i/j(V ; {Rj}) = A/, n(^o + Bf({Rj})), > 3 

i=l 



(2.3a) 



which resemble those for the su(2) generators Y a but with polynomial structural functions 
ip(Vo] {Rj}) depending additionally on invariant operators Ri : [V^,i?j] = unlike the 
quadratic function 

ip(Y ;J) = (Y + J) (J- y + l) (2-36) 

for su(2) [9]. Therefore, V a can be considered as generators of PLA su pd (2) acting on 
L{H) complementarily to Gi [5], i.e., [Gi, su pdi (2)] = and L(H) are decomposed into 
direct sums 

L(H)=^L([U])M[U]) = 8jm{\[l^^ 

U 

Ri\WJ) = k\UJ), N(f; [y) = / rj[^o + /-r)]- 1 / 2 , N(0; [/,]) = 1, VL|[/ i ]> = 0, 

r=0 

(2.4a) 

(Mf\[%f) = 8 [m 6 f ,,, / = E 7 h = E IM;/>«/l (2- 4fo ) 

m pi],/ 

of Gj ® SM p( ;(2)-irreducible subspaces J([/j]) with finite dimensions d([Zi]) where G 
L(H) are "lowest" weight basic vectors, quantum numbers l^i — 0, 1, . . . are eigenvalues 
of operators _Rj and "lowest" weight operators R are determined by the su pd {2) Casimir 
operators 

*(i? ) = ^(Vb; {#,}) - VV^_ = ^(Vb + 1; {Rj}) - V_V + , [V a , M>(i? )] = (2.4c) 

which, by construction (due to the invariant theory [5]), satisfy the characteristic identities 

*(Ro)\ HH ) = *(Z ) =^(l ;{lj}) = = ^(lo + ddh});^}) = *(Zo + d([ZJ)) (2.4d) 

which are useful in applications [5-9]. 

In the general cases the Gj-invariant Jordan mappings can be given as follows [9]: 

V+ = a+... a+a , VI = a h . . . a in a+, V = (E™ x Ni ~ N o)/( n + *)> 

R i = (EZi Ni + niV o)/( n + X )' ^i<i<m = H N i-i ~ Vj-i N 3 (2-5a) 
for models (2.1a) (/ij is the multiplicity of the a+ appearancein An ) and 

i M M 

V = J , V + = J + (aO™, VI = J_(a+r, J = -E <To(i), -4 = E ff ±«, 

1 i=i i=i 

i?i = Jo + (a+ai)/n, i? 2 = J, J(J + 1) = J 2 + (J+ J_ + J_ J+)/2 (2.56) 
for the model (2.1b). 

Appropriate structural polynomials (2.3a) are determined from Eqs. (2.5), (2.4c)- 
(2.4d). Specifically, they are given by expressions [9] 

il> D (V ; R u R 2 ) = (R 2 + V )(R 2 + 1 - VoXnR! +n- nV ) (n \ (2.6a) 



for the model (2.1b) and spread specifications of the model (2.1a) with H mps (l;n) = 
H hg (n) (high-harmonics-generation models) respectively; here (A)^ = A(A — 1) . . . (A — 
m+1). In the general case for the model (2.1a) we get fairly complex expressions (see, e.g., 
them for H mps (n;n) EE HH n) in [9]), and, therefore, a special (based on the Symbolic 
Computer Algebra System [14]) routine was developed to get such expressions automati- 
cally. 

The introduction of the collective operators above enables us to express all Hamiltonians 
(2.1) as follows [5,9]: 

H = h[AV + gV + + g*V- + C({^})], [V a , C] = 0, (2.7) 

where coefficients A and functions C = C({Ri}) are found from Eq. (2.1) with the help 
of Eqs. (2.5) for Ri, V a . In such a manner we get [9]: 

C hg (Ri) = (uji + cj )-Ri, A^ 9 = nuoi — coo for H h9 (n), (2.8a) 

Tl / Tb — 1 \ TL — 1 71 TL 

C^{[Ri\) = luoRi + E - " E iRi+i + E ^ E R ^ A/C = X> ~ "o, (2.86) 

3=1 71 \ i=l ) 3=1 i=j+l i=l 

C D (R 1 ,R 2 ) =nR 1 u 1 , A D = e-cui (2.8c) 

for H h9 (n), H^ c {n) = H mps (n; n) and H D respectively. At the same time decompositions 
(2.3) of L(H) are specified by determining quantum numbers k (integral of motions) and 
"lowest" weight vectors |[/j]) G L(H) as solutions of defining relations in (2.4a) [4,5]. 
Specifically, in such a manner one gets [9]: 

L h9 iU) =L«(lo = ^,h = ^) , IftOM)]) = ^^T/FlO), < K > n-l,s> 0, 

(2.9a) 



L fc ([h]) EE L fc 



l = (E K i ~~ S )l ( 1 + U )Jl = (E K i + nS )l ( X + n )> h<i<n = - «i 



i=l 



\M*i,8)]) = lli T' { * }y 10), «i>0, S >0, n«i = «. (2-96) 



k iV 

L D ([y) EE L D (l = -J = -l 2 ,h = -j + -;{lr>2+r=jr}), K > 0,j = — 

71/ Zi 



N' 
~2 



N 

' Y' 



\[k<2(K,j)};{jr}) = \K)f®\j;{3r};-j)a, |«)/ = M^WIO) G L F (1) (2.9c) 

where |j; {j r }', — j) a are "lowest" weight vectors of the s«(2)-irreducible "atomic" sub- 
spaces L(j;{j r }) C L a and extra integral of motions ("intermediate atomic quasispins" 
[3,9]) j r are due to the factor Sm in Gf\ in all these cases d([li\) = s + 1, where 
s = min{2j, n/n} for H D . 

So, Eqs. (2.4)-(2.9) yield the "s-u p d(2)-cluster" formulation of models (2.1) entailing a 
complete description of model dynamics with the help of the su p d(2) algebra formalism. 
In particular, it provides adequate representations [5-9] for model evolution operators 
Unit) to calculate quantum expectations (0(t)) of any operators O = 0(af ', . . .) acting 
onL(H): 

(0(t)) = Tr[U H (t)pV* H (t)0] = Y.Y.(UJ\pU\ I {t)OU H {t)\UJ). (2.10) 



3 Diagonal representations of model evolution oper- 
ators and algebraic schemes to determine them 

One of adequate representations of Unit) has a diagonal form [8,9] 

U H (t) = £ e ~ MUv) \ E v(lk})) (E v ([k})\, E v ([h}) = he([h];v) = h[C({h}) + K([h})} 

[k];v 

(3.1) 

where E v ([li\) and {\E V ( [h]))} are, respectively, eigenvalues and complete sets of orthonor- 
malized eigenvectors of the Hamiltonian (2.7): 

H\E f ([h})) = EfimEfdh})), (E f (5 fr , I = £ \E V ([U]))(E V ([U])\ (3.2) 

[k],v 

In the "linear" case, when = 2 in (2.3a), ip(Vo; {Rj}) = VK^oi J) an d PLA su p d{2) are 
reduced to the familiar su(2) algebra, the eigenproblem (3.2) is solved exactly [8] with 
the help of the SU(2) displacement operators = rg/\g\) = exp(£V + — £*F_) [11]; 

herewith solutions are given by simple analytical expressions [5-8]. 

However, it is not the case when > 3 in (2.3a) in view of the absence of explicit 
expressions for matrix elements ([/j]; f\ exp(J2i Q>iVi)\[li\; v) [8]. In fact, in [6] an algorithm 
has been developed to solve eigenproblem (3.2) with n$ > 3 via "dressing" an auxiliary 
"linear" (with = 2 one; but it is unsuitable for practical calculations, and, really 
nowadays it is known only an approximate analytical solution [8] of the problem (3.2) 
given by approximate SU (2)-quasiclassical eigenfunctions 

\Ef{[lU)) = \E^ 2 \Ui = rg/\g\)) = ^Y + - ^Y-)\[%v) =E^ B (OIW;/>, 

/ 

Sj v (0 = (9/\9\) f - v di j+f; _ j+v (2r), 2j = s = d[h] - 1 (3.3a) 

and eigenenergies 

EfiUi) = (E*([U]rf)\H\EF(MZ)) = *[C({li}) + K c (U,r)}, K c ([%r) = A(j + l ) 

-A(j - v) cos2r + 2\g\Y, j{s -/)(/ + l)<f>(-j + f) dP_ j+ft _ j+v (2r) <H j+f+h _ j+v (2r) 

f=o 

« Xt mf ([h];r) = A(j + Z ) - A(j - v) cos2r - 2\g\ (j - v) sin2r j<f>[{-j + v) cos2r], 

<K-j + f) = ([y ; /|V(n + i;R})M^o + i;i)IW;/) (3-36) 

where generators Y a of the su(2) algebra are connected with those of PLA su pd {2) 
via generalized Holstein-Primakoff mappings [5,8], d? mn {2r) is the SU{2)d- function ex- 
pressed in terms of the Gauss hypergeometric function [15] and approximate values 
-M5 (Mi r ) are calculated in the cluster mean-field approximation: ([/;]; /|F(V^)|[/j]; /) = 
F(([/j]; f\ V a \ /)). Values of the parameter r in (3.3b) are found from energy- st at ionarity- 
conditions and/or from minimizing a proximity measure of Hamiltonians H and if 9C (£) = 
J2 v ,[h] E$ c ([li\; £)\E$ c ([li\; £))(E% c ([li] ; £)|; a standard measure for such estimates on the 
subspaces L([li\) is defined with the help of the euclidean operator norm as follows [5,8] 

,2„,n T^jH-H^y _EA(K([k])) 2 -(K c (U,r)) 2 ] 



However, according to the general quasiclassicality theory [16] all approximations (3.3) 
are valid only for large values of d( [/«]), and, besides, the measure (3.4a) gives only a global 
rather than local characteristic of the approximate energy spectra {E^ c ([li\; £)} that does 
not allow to feel their important symmetry properties and local peculiarities related to 
"energy errors" 

AE v ([l i ])=h[X v ([k]) ~ \T(U,r)] = 5E v ([k])-E v ([k]) (3.46) 

Therefore, it is useful to have convenient algorithms to get exact solutions of the eigen- 
problem (3.2) using su p d(2) defining relations. 

Specifically, taking into account Eqs. (2.4a), one can look for eigenfunctions |P/([/j])) 
on each s-u p(i (2)-irreducible space £([/«]) in the form [5,8] 

\mU])) = E <#([»];«> = E Q f M\) V v + \[k\), Ql&i]) = N(v; M)Ql{[k\) (3.5a) 

v=0 v=0 

where s = d[li] — 1 and, in view of Eqs. (3.2) amplitudes Q{,([k]) satisfy the following 
orthonormalization and completeness conditions: 

E^TO Qv([h}) = 6 ff ,, E^(W) QUM) = (3.56) 

v=0 f 

Then, inserting Eq. (3.5a) for |Pj([/j])) and Eq. (2.7) for H in the first equation of (3.2) 
and using Eqs. (2.3a), (2.4c)-(2.4d), one gets a set of recurrence relations 

g *^l + v + 1; [k])Ql +1 (U) = [\f([k]) - A(v + l )]Q f v ([k})) - gQUdk}), v, f = o, • • • , s 

(3-6) 

which together with boundary conditions Q_\ ([/«]) = = Q s+ i([h]) determine amplitudes 
Ql([k]) and eigenenergies Ef([li\) from solutions of the Sturm-Liouville spectral problem 
[5,6] 

P V+1 (X) = [A - A(v + l )]P v (X) - \g\ 2 4>(l + v- [l^P^X), v = 0, • • • , s; (3.7a) 

P (A) = 1, P_i(A) = = [A - A(s + l )]Ps(X) - \g\ 2 ^(lo + s; [h})P s -i(X) (3.76) 
for finding non-classical orthogonal (in view of (3.5b)) polynomials 

P V (X) = (gTN- 2 (v; [k])Q v (U A)/Qo(M I A) = (gTN-^v; [k})Q v ([% A)/Qo(M I A) 

(3.8) 

for finding non-classical orthogonal (in view of Eqs. (3.5b) polynomials in the discrete 
variable A on the non-uniform lattice {A/([/j])}* =0 [5]. Indeed, Eqs. (3.7), (3.8) provide 
the following (easily realized by means of FORTRAN programs) algorithm for solving the 
eigenproblem (3.2). 

i) Using the recursive formula (3.7a) with initial values from Eq. (3.7b) one calculates 
the polynomial sequence {P„(A)}* =0 . 

ii) Inserting P S _ 1 (A), P S (A) in the last equality in (3.7b) one gets the algebraic equation 
with respect to A; its solution yield the sequence {A/([/j]) = A/}j =0 of admissible values 



iii) For each value A/([/j]) using {P V (X)}1 =0 and Eq. (3.8) one finds the sequence 
[Qli\h]) = Qv([k]] A/)}* =1 of all amplitudes as functions in the only undetermined quan- 
tity Q (\f) = Qo([k]',Xf) which, in turn, is found from the normalization condition of 
Eqs. (3.5b). 

In order to make numerical calculations we implemented this algorithm with the help 
of the REDUCE Package [14] similar to the algebraic construction in [17]. 



4 Numerical analysis for the second-harmonic gener- 
ation model 

In order to examine the efficiency of calculation schemes and the algorithm given above 
we tested them by means of computer experiments for the resonance second-harmonic- 
generation model widely examined in quantum optics [1,8,10] and determined by H hg (l; 2) 
with uj = 2uj\. In this case, according to Eqs. (2.6b), (2.8a), (2.9a), and (3.3b) we have 

V>(Zo + v + 1; [k]) = (k + 2v + 2)(k + 1 + 2v){s - v), k = 0, 1; s = 0, 1, . . . , (4.1a) 

C(h = !^-)=u 1 (k+2s), A = 0, XT f ([k];r) = -2\g\(j-v) sin2r y/flH + v) cos2r], 

4>{-j + f) = 2[s + 2k + 1 - (s - 2/)], j = s/2 (4.16) 

We implemented calculations of exact values Xf([k,s]) according to the algorithm 
of Section 3 and of their approximations X c ^^ {{k, s]; r^) according to Eq. (4.1b) for 
g = 1, k = 0, 1, s = 20, 100, 500, 1000, 10000. Values of the fitting parameter r were 
determined from energy-stationarity-conditions: 2r\ = arccos | (optimizing only the up- 
per part of spectra) [8], 2r 3 = arccos = | (quasi-linear approximation) [5,8] and from 
minimizing the proximity measure (3.4a): 2r 2 = arccos 4j ("smooth" cluster mean- field 

approximation) [8]; herewith X c p^ \[k, s]; =Fri) means that we take r = — n in the first half 
of spectra and r = r\ the second one. To estimate the accuracy of approximations we 
used, besides Eq. (3.4a) the measures 

X 2 m n_ ^=o[MM)-K mf (U,r)] 2 A2 nin _ E s v=s/2 [(X v ([h])-XT f ([h};r^ 
E[m ~ ZU(X v ([k})) 2 ' Eup[[ lD " EU /2 (A„(PJ)) 2 

(4.2) 

which characterize energy spectra more precisely in comparison with Eq. (3.4a). 
Typical results of numerical calculations obtained are presented in two tables. 

Table 1. Multiplets {A/([0, s])}, {A} m/ ([0, s}; n)} as functions in / = 0, 1, s for s = 10 2 



V 


A/(M) 




X c p f ([0,s]; T n) 


A^([0, S ];r 2 ) 


A^ / ([0, a ];r 3 ) 





-1536.9 


-1096.7 


-1545.3 


-1482.4 


-1421.2 


10 


-1151.7 


-919.6 


-1205.2 


-1175.2 


-1137.0 


20 


-798.1 


-720.0 


-880.0 


-873.3 


-852.7 


30 


-480.3 


-499.3 


-570.2 


-576.7 


-568.5 


40 


-205.5 


-259.0 


-276.7 


-285.6 


-284.2 



60 


205.5 


276.7 


276.7 


280.0 


284.2 


70 


480.3 


570.2 


570.2 


554.3 


568.5 


80 


798.1 


880.0 


880.0 


822.8 


852.7 


90 


1151.7 


1205.2 


1205.2 


1085.5 


1137.0 


100 


1536.9 


1545.3 


1545.3 


1342.3 


1421.2 


*H = 




10.222 


-12.220 


0.010 


-1.000 


*l = 




2.563 


0.670 


0.806 


0.657 


°E uP - 




0.670 


0.670 


0.944 


0.657 



Table 2. Multiplets {A/([0, s])}, {A} m/ ([0, s];^)} as functions in / = 0, s for s = 10 4 



v 


A/fro si) 




A c f m/ ([0 sl-^n) 


XT f (\0 si to) 


\T f (\0 sVr*) 





-1539573 


-1088743 


-1539658 


-1421266 


-1414284 


200 


-1460154 


-1055592 


-1470663 


-1364145 


-1357712 


600 


-1304855 


-986414 


-1334422 


-1249971 


-1244570 


1000 


-1154419 


-913501 


-1200541 


-1135886 


-1131427 


2000 


-800913 


-715586 


-876392 


-851069 


-848570 


3000 


-483083 


-496534 


-567801 


-566815 


-565713 


4000 


-208106 


-257638 


-275425 


-283125 


-282856 


4400 


-113484 


-156775 


-163176 


-169807 


-169714 


4800 


-31728 


-52979 


-53690 


-56579 


-56571 


5000 

















5200 


31728 


53690 


53690 


56557 


56571 


5600 


113484 


163176 


163176 


169603 


169714 


6000 


208106 


275425 


275425 


282559 


282856 


7000 


483083 


567801 


567801 


564553 


565713 


8000 


800913 


876392 


876392 


845978 


848570 


9000 


1154419 


1200541 


1200541 


1126836 


1131427 


9400 


1304855 


1334422 


1334422 


1239020 


1244570 


9800 


1460154 


1470663 


1470663 


1351113 


1357712 


10000 


1539573 


1539658 


1539658 


1407125 


1414284 


*Ir = 




11.102 


-11.122 


0.000 


-0.010 


*! = 




2.577 


0.590 


0.626 


0.625 


/> 2 

°E u p - 




0.590 


0.590 


0.639 


0.625 



As is seen from data given in Tables 1, 2, there is a satisfactory consent of exact results 
and those obtained with the help of the approximate formula (3.3b) at s ^> 1 in the most 
parts of energy spectra for r = rj =2 3 and r = =pri. Discrepancies between exact and 
approximate results in middle parts of spectra are explained by the availability of the 
square-root singularities in H cm f(0 = £,,[y E^([k}; £)\E™f {[% 0)(K mf &i] 5 01, and 
the symmetry breaking of spectra at r = r i= i i2 is due to smoothing this singularities 
within the framework of the "smooth" su(2) - quasiclassical approximation [8,9]. A more 
full analysis of results of computer experiments above will be given elsewhere. 



5 Conclusion 



So, we developed an algebraic scheme and routines (implemented with the help of the 
REDUCE Package [14]) to get exact solutions of both spectral and evolution problems 
for a wide class of integrable models of quantum optics having PL A su p d{2) as their 
dynamic symmetry algebras. These results along with appropriate computer realizations 
of Eqs. (3.3) provide an efficient tool for examining models under study in appropriate 
(determined by computer powers) ranges of characteristic model parameters given by Zj. 
Numerical calculations given in Section 4 showed a good consent of exact and approximate 
results at s ^> 1 and at relevant choices of the fitting parameter r in (3.3b). Furthermore, 
the approximate solutions of the eigenproblem (3.2) can be improved with the help of 
the algebraic perturbation algorithm like that developed in [17] for the hydrogen atom 
or by means of using modifications of algebraic schemes of the work [6]. Therefore, 
the exact calculation schemes developed above, when being completed by quasiclassical 
calculations [8] based on formulas like Eq. (3.3) as well as their improvements, can be 
used for analyzing models under consideration in all ranges of characteristic parameters 
/j. The work along these lines is in progress. References 
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